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Production of perfect crystals, free of residual strain and dislocations and 
with prescribed dopant concentration, by the Czochralski method is possible 
only if the complex, interacting phenomena that affect crystal growth in a 
Cz-puller are fully understood and quantified. Natural and forced convection 
in the melt, thermocapillary effect and heat transfer in and around the 
crystal affect its growth rate, the shape of the crystal-melt interface and 
the temperature gradients in the crystal. In this work we have 
concentrated on describing the heat transfer problem in the crystal and 
between the crystal and all other surfaces present in the crystal pulling 
apparatus . 

This model and computer algorithm are based on the following assumptions: i) 
only conduction occurs in the crystal (experimentally determined conductivity 
as a function of temperature is used), ii) melt temperature and the 
melt-crystal heat transfer coefficient are available (either as constant 
values or functions of radial position), iii) pseudo-steady state is achieved 
with respect to temperature gradients, iv) crystal radius is fixed, v) both 
direct and reflected radiation exchange occurs among all surfaces at various 
temperatures in the crystal puller enclosure. 

This is the first model to replace the simple Stefan’s law for radiation from 
the melt and crystal surface to the ambient with the Gebhart’s enclosure 
theory which accounts for both direct and reflected radiation. Available 
emissivities were used. The position of the melt-crystal interface is not 
assumed a priori but is computed as part of the solution as a locus of points 
where the flux continuity (heat flux from the melt plus the rate of release of 
the heat of crystallization is balanced by heat conduction into the crystal) 
is satisfied at the melting point isotherm. 

Galerkin finite elements with triangular basic cells and linear trial 
functions are used to reduce the problem to a set of nonlinear algebraic 
equations. These are solved by an iterative scheme both for the temperatures 
and the melt-crystal interface position. 

Using 10 elements in the radial direction and 30 in the axial proved adequate 
to achieve the desired accuracy and resulted in 352 unknowns. Doubling the 
number of nodes did not affect the results. One minute of CPU time on a DEC20 
computer was required per iteration. Convergence was obtained usually in less 
than 8 iterations. 

The model calculates the axial and radial temperature profiles in the crystal 
and the position of the melt-crystal interface as a function of the input 
parameters (melt temperature, crucible wall temperature, temperatures of 
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various surfaces in the enclosure, pulling rate, crystal radius). Comparison 
with the results in the literature^) which used Stefan's law to describe 
radiative heat transfer show that these underpredict the axial temperature 
profile and the radial temperature gradients. Thus, proper accounting for 
reflected radiation seems important. The model also indicates that an 
increase in pulling rate at other fixed parameters would tend to make the 
melt-crystal interface more concave to the crystal (i.e. melt protrudes into 
the crystal). The same trend is caused by an increase in melt or crucible 
temperature . 

These final predictions regarding the melt-crystal interface need further 
verification. This can only be done by incorporating the thermocapillary 
effect, to describe the meniscus at the crystal-melt-gas line of contact, and 
by evaluating the convective patterns in the melt to obtain proper melt 
temperature and heat transfer coefficient profiles. This work is currently in 
progress . 


d) Kobayashu, N., Heat Transfer in Cz Crystal Growth (W.R. Wilcox, ed.), 
Marcel Dekker, Inc. (1981). 

*2) Ramachandran, P.A. and Dudukovic, M.P., "Simulation of Temperature 
Distribution in Crystals Grown by Czochralski Method," J. Crystal Growth 
(accepted for publication, 1984). 

(3) Williams, G. and Rousser, P.E., J. Crystal Growth, 64 , 448 (1983). 
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Fig. 3 — Schematic drawing illustrating the spatial relationship of 
separate factors that affect oxygen content in silicon during 
crystal growth. (After a drawing in Benson, et al. [25].) 


R. R. Huff, Solid State Technology, Feb. *83 
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Problem Complexities 


1. Temperature distribution in melt and crystal 
are coupled. 

2. Complex flow patterns in melt 

1. forced convection 

2. natural convection 

3. thermocapillary effect 

3. Floating boundary, interface shape cannot be 
fixed a priori 

4. Radiation interactions 
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Variables Affecting Cz Growth 


i) geometric variables including crucible shape, 

enclosure shape, heater position and shape, etc.; 

ii) controllable operating variables such as pulling 
rate, crystal rotation rata, crucible rotation 
rate, gas flow rate, furnace power; 

iii) process variables such as melt depth, crystal 
diameter, etc. 

Due to the complex interaction among many variables 
their interrelationship cannot be understood in a 
quantitative sense through any number of experiments 
unless experimental results are interpreted through an 
appropriate model of the system. 
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Ultimate Objective Is to Produce a Perfect Crystal 


CRYSTAL 


. FREE OF RESIDUAL STRAIN 
. FREE OF DISLOCATIONS 
. UNIFORMLY DOPED 

ONCE A PROPER MODEL IS ESTABLISHED AND 
VERIFIED THE EFFECT OF CHANGES IN VARIOUS VARIABLES 
ON CRYSTAL QUALITY CAN BE ASSESSED MORE READILY 
WITH FAR FEWER EXPERIMENTS. OPTIMIZATION OF THE 
OPERATION BECOMES POSSIBLE. 
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Crucible Heater 


Figure 1. Schematic of Czochralski single crystal pulling apparatus, 

indicates radiation interaction between various surfaces 



Enclosures for Crystal Problem 



CRYSTAL IS DIVIDED INTO NB2 ELEMENTS. THESE PROVIDE 
NB2 SURFACES. THREE ADDITIONAL SURFACES ARE CONSIDERED 
AS SHOWN ABOVE 


N = NB2 + 3 
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Model Equations 
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Calculation of Radiation Flux 

GEBHAR.T METHOD IS USED 


CONSIDER AN ENCLOSURE WITH N SURFACES, EACH SURFACE 
BEING AT TEMPERATURE T k . THE HEAT LOSS FROM SURFACE 
k IS 


Q k = A k € k c ' T 
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CALCULATION G| k (J = 1 TO N) FOR ANY GIVEN SURFACE k 




SURFACE k RECEIVES BOTH DIRECT AS HELL AS REFLECTED 
RADIATION FROM SURFACE j. DIFFUSE REFLECTION IS 
ASSUMED. 
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DIRECT RADIATION (j to k) 
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REFLECTED RADIATION (j to n to k) 
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F|j - 0 for I < NB2 
and 

J < NB2 

HENCE THE CALCULATION OF 6 jk CAN BE SIMPLIFIED. 
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LINEARIZED FORM OF BOUNDARY CONDITION 
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AN EFFECTIVE TEMPERATURE FOR RADIATION CAN BE DEFINED 
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Finite-Element Equations 



and P depend on nodal coordinates, physical properties., 
and the boundary conditions imposed on the element. 

Equations for each element can be assembled and the 
resulting global matrix can be solved for all the nodal 
temperatures. 
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TO START ITERATION SCHEME USE ONE -DIMENSIONAL MODEL 
WITH FINITE PECLET NUMBER 

(Wilcox & Duty, J. Heat Transfer, 1966) 




PULLING RATE EFFECT NEGLIGIBLE 
(SMALL PECLET NUMBER) WHEN 



8hR 
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Program Structure 


1. INPUT 

2. NODE AND ELEMENT NUMBERING 

3. ASSUME TRIAL VALUES FOR TEMPERATURES FOR FIRST 
ITERATION 

3a. CALCULATE G FACTORS 
3b. CALCULATE 

T eff,R' h r' V etc * 

4. SOLVE FINITE ELEMENT EQUATIONS FOR. NODAL TEMPERATURES 

5. CHECK FOR CONVERGENCE 

6* T New = ^1 ^calculated 

+ (1 - W^) T assurneC j 

7. LOCATE Z POSITIONS WHERE TEMPERATURE IS EQUAL TO 
TO MELTING POINT (INTERFACE SHAPE). 

8. REPEAT CALCULATIONS TILL CONVERGENCE IS OBTAINED. 
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Special Features 


. INTERFACE SHAPE CALCULATED AND NOT ASSUMED 

. TWO DIMENSIONAL MODEL 

. ENCLOSURE THEORY PERMITS DETAILED ANALYSIS OF 
RADIATION HEAT TRANSFER 

. EFFECTS OF ADDITIONAL HEATING AND COOLING SURFACES 
CAN EASILY BE SIMULATED 

. PHYSICAL PROPERTIES ASSUMED AS A FUNCTION OF 
TEMPERATURE 
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Figure 2. Interface shapes c< 

geometric parameters ai 
element meshes used 





Table 2. Parameter Values Used in Illustrative Application Example 


radius of crystal 

height of crystal 

crystal pulling rate 

crucible radius 

height of exposed crucible 

wall area 

wall temperature 

melt temperature 

crucible temperature 

average argon temperature 

effective temperature for radiation 
(Stefans model) 

convective heat transfer coefficient 

melt to crystal heat transfer 
coefficient 

thermal conductivity 

emissivity of silicon crystal 

emissivity of melt 
emissivity of crucible 
emissivity of wall 


4 cm 
17 cm 
0.002 cm/s 
10 cm 
3 cm 
1500 cm 2 
100°C 
1470°C 
1480°C 
250°C 

327°C 

1.2433 x 10“ 4 (AT) 1/3 W/cm 2 K 
0.3265 W/cm 2 K 

0.98892408 - 9.4286595 x 10 -4 T 
+ 2.889 x 10“ 7 T 2 , W/cm K 
(T in Kelvin) 

£ * 0.64 if T < 1000 K 
= .9016 - 2.616 x 10" 4 T 
if T > 1000 K 
0.3 
0.59 
0.59 
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SURFACE TEMPERATURE 



AXIAL DISTANCE FROM MELT SURFACE, cm 


Figure 5. Predicted crystal surface temperatures based on Gebhart and Stefan models and experimental 
results of Williams and Reusser [3] 






AXIAL DISTANCE 




RADIAL DISTANCE, cm 


Figure 6 . Interface shape calculated for the data 
of Table 2. — Gebhart analysis for 
radiation, ■» — Stefan model 
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Summary 


A FINITE ELEMENT ALGORITHM HAS BEEN DEVELOPED 
FOR EVALUATION OF THE TEMPERATURE PROFILES IN THE 
CRYSTAL AND FOR PREDICTION OF THE POSITION OF THE 
CRYSTAL-MELT INTERFACE. 

IT IS SHOWN THAT ACCOUNTING PROPERLY FOR 
RADIATION EFFECTS CAN INFLUENCE CONSIDERABLY THE 
PREDICTIONS FOR THE TEMPERATURE PROFILE AND MELT- 
CRYSTAL INTERFACE, 
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SUMMARY (Cont'd) 


THE DEVELOPED ALGORITHM, AT PRESENT, INCORPORATES 
THE FOLLOWING SIMPLIFYING ASSUMPTIONS 

- pseudo-steady state 

- symmetry 

- constant melt temperature 

- constant heat transfer coefficient at 
melt-crystal interface 

- no thermo capillary effect 

- fixed crystal radius 

COUPLING OF THIS MODEL WITH A REALISTIC HYDRO- 
DYNAMIC AND THERMOCAPILLARY MODEL IS NECESSARY AND IS 
CURRENTLY IN PROGRESS, 

COMPARISON WITH EXPERIMENTAL DATA AND PARAMETER 
SENSITIVITY STUDIES ARE ALSO ESSENTIAL. 
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Recommendations 


OUR UNDERSTANDING OF THE COMPLEXITIES OF THE 
CZ (OR LECZ) PROCESS IS FAR FROM COMPLETE. 

UNIVERSITY PERSONNEL AT A NUMBER OF LOCATIONS 
MUST BE ENCOURAGED (AND FUNDED) TO CONTINUE THE WORK 
IN THIS AREA. 

COMPETITION AND OVERLAP OF COVERAGE OF CERTAIN 
TOPICS IS HELPFUL. INADEQUATE COVERAGE AND NEGLECT 
OF THESE PROBLEMS CAN ONLY HURT OUR COMPETITIVENESS 
IN THE FUTURE. 
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DISCUSSION 


WITT: I would like you to treat the inverse problem. Namely, we need an 

establishment of the boundary conditions, or the range of boundary 
conditions, under which in the desired thermal configuration the 
crystal in the melt can be achieved, rather than quantitative 
verification through modeling of adverse thermal configurations. I 
recognize that the inverse problem is infinitely more difficult to 
solve. What I think industry needs to know is, what are the 
required boundary conditions necessary to achieve optimized silicon 
growth? Not verification of the complexity of the heat transfer 
conditions through modeling. 

DUDUKOVIC: Your point is well taken, but I think the two really go together. 

1 don't see it at all incompatible in the process of learning the 
complexities, that you show that one should be thinking about some 
other configurations in which asymmetry will not be a problem. Even 
if you go to those other configurations, it's still nice to have a 
quantitative rigorous model for well-defined situations. 

WITT: I'm only trying to preserve realism, and the realism is that you can 

only sell the effort to industry if industry sees the direct 
impact. I think industry wants concrete input to systems design and 
the connection from what we are doing to that is obscure. 

SCHWUTTKE: It looks like you are trying to calculate the position of the train 

after it left the station. The crystals are already grown and 
nothing happens anymore that I am concerned about. I need to know 
how to shape my interface. 

DUDUKOVIC: In order to know where the interface shape is — you are not going 

to get there just by knowing the meniscus. We also need to know the 
radiative heat fluxes and the conductive heat fluxes. 

One has to start from someplace. Once you have solved the crystal 
problem right you can now couple that with any level of 
approximation for the meniscus and for the melt that you want. The 
three are coupled. You can not arrive at a one- line formula by 
which you can resolve the interface shape just by considering one of 
those effects. 

MORRISON: I know for a fact that industry has been interested in your 

modeling. I think that when the models are built and verified you 
can go back and use the models as a tool to build a station. 

Have you compared the interfaces of the terminations of grown 
crystals with the interface shapes that your model is generating? 

DUDUKOVIC: That would be the wrong thing to do, because it would just 

demonstrate the effect of one factor on the position of that 
interface. This is exactly the trap that industrial people often 
fall into. They run out and say it doesn't match my experimental 
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observations. One can make sense of complex phenomena by putting 
them together and the only way that putting them together can be 
done is through mathematics. There is no other way. Everything 
else is just handwaving. 
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